Distribution Map: Bombus and its Subgenera

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

library('rgdal')
library('dismo')
library('raster')
library('sp')
library('maptools')
library('readr')
library('ggplot2')
library('ggmap')
library('sf')

#saving file into a label and reading
occ.data <- read.csv(file = "data/bombus_emu.csv")

summary(occ.data)

            Barcode     
 USNMENT00123087:    1  
 USNMENT00679000:    1  
 USNMENT00679001:    1  
 USNMENT00679002:    1  
 USNMENT00679003:    1  
 USNMENT00679004:    1  
 (Other)        :44108  
                                                                                        Taxon      
 Bombus unidentified : Apidae : Hymenoptera : Insecta : Arthropoda                         :10752  
 Bombus bifarius nearcticus : Apidae : Hymenoptera : Insecta : Arthropoda                  : 2533  
 Bombus (Pyrobombus) mixtus Cresson, 1878 : Apidae : Hymenoptera : Insecta : Arthropoda    : 2122  
 Bombus (Bombus) affinis Cresson : Apidae : Hymenoptera : Insecta : Arthropoda             : 1902  
 Bombus (Pyrobombus) flavifrons Cresson, 1863 : Apidae : Hymenoptera : Insecta : Arthropoda: 1701  
 Bombus (Pyrobombus) bifarius Cresson, 1878 : Apidae : Hymenoptera : Insecta : Arthropoda  : 1553  
 (Other)                                                                                   :23551  
 DarDecimalLatitude DarDecimalLongitude
 Min.   :-53.63     Min.   :-172.72    
 1st Qu.: 39.52     1st Qu.:-122.42    
 Median : 42.08     Median :-105.27    
 Mean   : 43.63     Mean   : -94.24    
 3rd Qu.: 48.09     3rd Qu.: -84.19    
 Max.   : 86.33     Max.   : 175.38    
 NA's   :18110      NA's   :18110      
                      DarScientificName DarScientificNameAuthor
 Bombus unidentified           :10752             :13285       
 Bombus bifarius nearcticus    : 2533   Cresson   :11766       
 Bombus (Pyrobombus) mixtus    : 2122   Smith     : 3193       
 Bombus (Bombus) affinis       : 1902   Nylander  : 2079       
 Bombus (Pyrobombus) flavifrons: 1701   (Linnaeus): 1868       
 Bombus (Pyrobombus) bifarius  : 1553   Greene    : 1781       
 (Other)                       :23551   (Other)   :10142       
 DarScientificNameAuthorYear    DarSpecificEpithet
              :13285         unidentified:10752   
 Cresson, 1878: 4436         bifarius    : 4086   
 Cresson, 1863: 4097         mixtus      : 2122   
 Smith        : 3080         affinis     : 1902   
 Cresson      : 2861         flavifrons  : 1701   
 Nylander     : 2079         melanopygus : 1423   
 (Other)      :14276         (Other)     :22128   
 DarYearCollected       DarCollector  
 Min.   : 947     L. Macior   :15123  
 1st Qu.:1932     [Not Stated]: 4316  
 Median :1965     D. Graham   : 4200  
 Mean   :1953     J. Bridwell : 2096  
 3rd Qu.:1971     E. Stiles   : 1582  
 Max.   :2017     J. Clarke   : 1523  
 NA's   :4567     (Other)     :15274  
                                      AdmGUIDPreferredValue
 ark:/65665/300030cb7-f6da-4dc3-8005-b769a0386842:    1    
 ark:/65665/30003be57-6cf7-4e31-b1e6-df2a59f21999:    1    
 ark:/65665/3000427d0-37c4-4d9e-8e40-494c334429df:    1    
 ark:/65665/30004595f-cd68-4adb-9cb1-50905d58f9e0:    1    
 ark:/65665/30004b971-c32d-44a6-b1e1-c774d8f0da37:    1    
 ark:/65665/3000513bb-73ae-48ed-bfe3-f47fa9a4cf02:    1    
 (Other)                                         :44108    
        DarStateProvince        DarCounty    
 Colorado       : 5326   [Not Stated]:15382  
 Yukon Territory: 3953   Boulder     : 3525  
 Washington     : 3496   Tuolumne    : 2025  
 Sichuan        : 3428   Racine      : 1169  
 California     : 3378   Warren      : 1139  
 [Not Stated]   : 2834   Surrey      :  811  
 (Other)        :21699   (Other)     :20063  
       DarLocality    DarMonthCollected DarDayCollected
 [Not Stated]:14646   Min.   : 1.000    Min.   : 1.0   
 Kluane      : 3771   1st Qu.: 6.000    1st Qu.: 8.0   
 Esher       : 1080   Median : 7.000    Median :15.0   
 Clearview   :  472   Mean   : 6.898    Mean   :15.5   
 Corvallis   :  434   3rd Qu.: 8.000    3rd Qu.:23.0   
 Mount Emei  :  415   Max.   :12.000    Max.   :31.0   
 (Other)     :23296   NA's   :4910      NA's   :7562   
 DarFieldNumber DarCollectorNumber DarPreparations
 Mode:logical   Mode:logical       Pinned:44114   
 NA's:44114     NA's:44114                        
                                                  
                                                  
                                                  
                                                  
                                                  

dir.create(path = 'data')

'data' already exists

dir.create(path = 'output')

'output' already exists

#shows properties of data, rows and columns
dim(occ.data)

[1] 44114    19

#Column names
colnames(occ.data)

 [1] "Barcode"                     "Taxon"                      
 [3] "DarDecimalLatitude"          "DarDecimalLongitude"        
 [5] "DarScientificName"           "DarScientificNameAuthor"    
 [7] "DarScientificNameAuthorYear" "DarSpecificEpithet"         
 [9] "DarYearCollected"            "DarCollector"               
[11] "AdmGUIDPreferredValue"       "DarStateProvince"           
[13] "DarCounty"                   "DarLocality"                
[15] "DarMonthCollected"           "DarDayCollected"            
[17] "DarFieldNumber"              "DarCollectorNumber"         
[19] "DarPreparations"            

#bioclimatic data
bioclim.data<- getData(name = 'worldclim',
                       var = 'bio',
                       res = 2.5,
                       path = "data/")

#disregarding all nulled values for lon and lat
clean.data <- subset(occ.data, !is.na(DarDecimalLatitude) & !is.na(DarDecimalLongitude))
#shows properties of data

dim(clean.data)

[1] 26004    19

clean.data

NA

#Determining what the max and min geographical values of map are
max.lat <- (max(clean.data$DarDecimalLatitude))
min.lat <- (min(clean.data$DarDecimalLatitude))
max.lon <- (max(clean.data$DarDecimalLongitude))
min.lon <- (min(clean.data$DarDecimalLongitude))

#Display Values
max.lat

[1] 86.3264

min.lon

[1] -172.72

max.lon

[1] 175.38

min.lat

[1] -53.63

#Establishing extent of map (will be able to plot only north american in later model)
geo.extent <- extent(x = c(min.lon, max.lon, min.lat, max.lat))
c(min.lon, max.lon, min.lat, max.lat)

#Using only collumns from lon and lat
loc.data <- clean.data [, c('DarDecimalLongitude', 'DarDecimalLatitude')]
loc.data

NA

data("wrld_simpl")

plot(wrld_simpl,
     xlim = c(min.lon, max.lon),
     ylim = c(min.lat, max.lat ), 
     axes = TRUE,
     col = 'grey95')

#adding points on to map
points(x = loc.data$DarDecimalLongitude,
       y = loc.data$DarDecimalLatitude,
       col = 'olivedrab',
       pch = 20,
       cex = 0.75)

#box around map
box()

bc.model <- bioclim(x = bioclim.data, p = loc.data)

predict.presence <- dismo::predict(object = bc.model, x = bioclim.data, ext = geo.extent)

#plot base map
plot(wrld_simpl,
     xlim = c(min.lon, max.lon),
     ylim = c(min.lat, max.lat),
     axes = TRUE,
     col = "grey95")

#adding model probabilities
plot(predict.presence, add = TRUE)

 
#redrawing country borders
plot(wrld_simpl, add = TRUE, border = "grey6")

#adding OG observations
points(x = loc.data$DarDecimalLongitude,
       y = loc.data$DarDecimalLatitude,
       col = 'olivedrab',
       pch = 20,
       cex = 0.75)

NA
NA

sep.data <- clean.data [, c('Taxon', 'DarDecimalLongitude', 'DarDecimalLatitude')]

pyro.data <- read.csv(file = 'data/Pyrobombus_emu.csv')

locpyro.data <- pyro.data [, c('DarDecimalLongitude', 'DarDecimalLatitude')]

locpyro.data

plot(wrld_simpl,
     xlim = c(min.lon, max.lon),
     ylim = c(min.lat, max.lat),
     axes = TRUE,
     col = "grey95")

plot(wrld_simpl, add = TRUE, border = "grey5")

points(x = clean.data$DarDecimalLongitude,
       y = clean.data$DarDecimalLatitude,
       col = 'olivedrab',
       pch = 20,
       cex = 0.75)

points(x = locpyro.data$DarDecimalLongitude,
       y = locpyro.data$DarDecimalLatitude,
       col = 'red',
       pch = 20,
       cex = 0.75)

Agro.data <- read.csv(file = 'data/Agro_emu.csv')
Alpigeno.data <- read.csv(file = 'data/Alpigeno_emu.csv')
Alpino.data <- read.csv(file = 'data/Alpino_emu.csv')
Bombias.data <- read.csv(file = 'data/Bombias_emu.csv')
BombusSub.data <- read.csv(file = 'data/BombusSub_emu.csv')
Brachy.data <- read.csv (file = 'data/Brachycephalobombus_emu.csv')
Coccineo.data <- read.csv(file = 'data/Coccineo_emu.csv')
Confusii.data <- read.csv(file = 'data/Confusii_emu.csv')
Crotchii.data <- read.csv(file = 'data/Crotchii_emu.csv')
Collumano.data <- read.csv(file = 'data/Collumano_emu.csv')
Diverso.data <- read.csv(file = 'data/Diverso_emu.csv')
Dasy.data <- read.csv(file = 'data/Dasy_emu.csv')
Fervido.data <- read.csv(file = 'data/Fervido_emu.csv')
Thoraco.data <- read.csv(file = 'data/Thoraco_emu.csv')

#Disregarding all na values in each data frame
Alpigeno.data <- subset(Alpigeno.data, !is.na(DarLatitude) & !is.na(DarLongitude))
locAlpigen.data <- Alpigeno.data [, c('DarLongitude', 'DarLatitude')]
locAlpigen.data

Agro.data <- subset(Agro.data, !is.na(DarLatitude) & !is.na(DarLongitude))
locAgro.data <- Agro.data [, c('DarLongitude', 'DarLatitude')]

Alpino.data <- subset(Alpino.data, !is.na(DarLatitude) & !is.na(DarLongitude))
locAlpino.data <- Alpino.data [, c('DarLongitude', 'DarLatitude')]

Bombias.data <- subset(Bombias.data, !is.na(DarLatitude) & !is.na(DarLongitude))
locBombias.data <- Bombias.data [, c('DarLongitude', 'DarLatitude')]


BombusSub.data <- subset(BombusSub.data, !is.na(DarLatitude) & !is.na(DarLongitude))
locBombus.data <- BombusSub.data [, c('DarLongitude', 'DarLatitude')]


Brachy.data <- subset(Brachy.data, !is.na(DarLatitude) & !is.na(DarLongitude))
locBrachy.data <- Brachy.data [, c('DarLongitude', 'DarLatitude')]

#skipped Coccineo

Collum.data <- subset(Collumano.data, !is.na(DarLatitude) & !is.na(DarLongitude))
locCollum.data <- Collum.data [, c('DarLongitude', 'DarLatitude')]

Thoraco.data <- subset(Thoraco.data, !is.na(DarLatitude) & !is.na(DarLongitude))
locThoraco.data <- Thoraco.data [, c('DarLongitude', 'DarLatitude')]

plot(wrld_simpl,
     xlim = c(min.lon, max.lon),
     ylim = c(min.lat, max.lat),
     axes = TRUE,
     col = "grey95")

plot(wrld_simpl, add = TRUE, border = "grey5")

#points(x = clean.data$DarDecimalLongitude,
 #      y = clean.data$DarDecimalLatitude,
  #     col = 'olivedrab',
   #    pch = 20,
    #   cex = 0.75)

points(x = locpyro.data$DarDecimalLongitude,
       y = locpyro.data$DarDecimalLatitude,
       col = 'red',
       pch = 20,
       cex = 0.75)

points(x = locAlpigen.data$DarLongitude,
       y = locAlpigen.data$DarLatitude,
       col = 'yellow',
       pch = 20,
       cex = 0.75)

points(x = locAgro.data$DarLongitude,
       y = locAgro.data$DarLatitude,
       col = 'pink',
       pch = 20,
       cex = 0.75)
 
points(x = locAlpino.data$DarLongitude,
       y = locAlpino.data$DarLatitude,
       col = 'blue',
       pch = 20,
       cex = 0.75)

points(x = locBombias.data$DarLongitude,
       y = locBombias.data$DarLatitude,
       col = 'purple',
       pch = 20,
       cex = 0.75)

points(x = locBombus.data$DarLongitude,
       y = locBombus.data$DarLatitude,
       col = 'grey',
       pch = 20,
       cex = 0.75)

points(x = locBrachy.data$DarLongitude,
       y = locBrachy.data$DarLatitude,
       col = 'olivedrab',
       pch = 20,
       cex = 0.75)

points(x = locCollum.data$DarLongitude,
       y = locCollum.data$DarLatitude,
       col = 'black',
       pch = 20,
       cex = 0.75)

points(x = locThoraco.data$DarLongitude,
       y = locThoraco.data$DarLatitude,
       col = 'green',
       pch = 20,
       cex = 0.75)

plot(wrld_simpl,
     xlim = c(-170, -39),
     ylim = c(11, 72),
     axes = TRUE,
     col = "grey95")
points(x = locpyro.data$DarDecimalLongitude,
       y = locpyro.data$DarDecimalLatitude,
       col = 'red',
       pch = 20,
       cex = 0.75)

points(x = locBombus.data$DarLongitude,
       y = locBombus.data$DarLatitude,
       col = 'grey',
       pch = 20,
       cex = 0.75)
points(x = locCollum.data$DarLongitude,
       y = locCollum.data$DarLatitude,
       col = 'black',
       pch = 20,
       cex = 0.75)

points(x = locThoraco.data$DarLongitude,
       y = locThoraco.data$DarLatitude,
       col = 'green',
       pch = 20,
       cex = 0.75) 

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.